home *** CD-ROM | disk | FTP | other *** search
/ Info-Mac 4 / Info_Mac IV CD-ROM (Pacific HiTech Inc.)(August 1994).iso / Science / MAL ©P.f.Howden 1⁄1⁄89 / DIFFSIMP < prev    next >
Text File  |  1989-03-22  |  5KB  |  82 lines
  1. 1 CLS:CLEAR:PRINT"DIFFSIMP-Solves Groups of Initial value,Non/linear,SIMultaneous,IMPlicit,":PRINT"ordinary Differential Equations to any order,including related or unrelated ":PRINT"Algebraic equations."
  2. 2 PRINT"Test eqns in listing DIFFSIMPEQNS.TXT can be MERGED & appear below,OR"
  3. 3 PRINT"TYPE EQUATIONS ON LINES 10-19  at end of listing":PRINT"eg:  10 F=X(1,1)*X(0,2)+T:RETURN   ..........":PRINT"Press SPACE bar if iterations not converging."
  4. 4 PRINT"DEFINE:VARIABLES X(I,J)=X[identify number(>=0), differential order(>=0)].":PRINT"IF USING VARIABLES LIMITER,eg TYPE:  10 GOSUB 740:F=X(0,1)-T:RETURN":PRINT"But see text & REM on line 8. For all equations, F=0.
  5. 5 DEFDBL C-D,F-G,Q,V-Z:'--------------INPUTS+INITIAL CONDITS to line 250
  6. 6 PRINT:INPUT "NUMBER OF EQUATIONS=";N:INPUT "MAXIMUM ORDER=";E:INPUT"R (usually 1,more if difficult to converge)=";R:N=N-1:INPUT "ACCCURACY,AC (1E-08?)=";AC:ST=500:G=.25:H=.1
  7. 7 B=1:INPUT"STARTING To (typically=0) =";T0
  8. 8 'Equations are solved at T0+H/2.NOTE line 500 for  optional ending conditions.Press P if any P large,R not changing,X iterations barely changing.
  9. 9 GOTO 30
  10. 30 DIM S(N),S1(N),X(N,E),Y(N),P(N),R(N),C(N,E),O(N),A(N)
  11. 120 FOR J=0 TO N:IF E=0 OR N=0 THEN O(J)=E
  12. 123 IF E=0 THEN 165
  13. 127 IF N=0 THEN 140
  14. 130 PRINT"ORDER OF VARIABLE";J;:INPUT"=";O(J):IF O(J)=0 THEN 165
  15. 140 FOR I=0 TO O(J)-1
  16. 150 PRINT"INPUT INITIAL C(";J;",";I;")=";:INPUT C(J,I)
  17. 160 NEXT I:GOTO 170
  18. 165 PRINT"STARTING X";J;"=";:INPUT C(J,0)
  19. 170 A(J)=1+J:S(J)=-1:NEXT J
  20. 180 INPUT "ANY EQUATION CHANGEABLE CONSTANT M=";M
  21. 190 IF N=0 THEN 230
  22. 200 PRINT"SAME EQUATION SEQUENCE [currently";:FOR I=0 TO N:PRINT A(I)-1;:NEXT I:INPUT "], (1/0)";W:IF W=1 THEN 230
  23. 210 PRINT"EQUATION SEQUENCE:":FOR I=0 TO N:PRINT"EQUATION";I;"->";:INPUT A(I)
  24. 220 A(I)=A(I)+1:NEXT I
  25. 230 PRINT"SAME SIGN/S  [currently";:FOR I=0 TO N:PRINT S(I);:NEXT I:INPUT "], (1/0)";W:IF W=1 THEN 250
  26. 240 FOR I=0 TO N:PRINT"SIGN (±1) FOR EQUATION";I;"=";:INPUT S(I):NEXT I
  27. 250 PRINT"INPUT T INCREMENT  H(currently";H;")=";:INPUT H
  28. 260 '------------------------------------PROGRESS DECISIONS to line 350
  29. 270 IF U=0 THEN INPUT "MENU: 0=AUTO: 1=INCREMENT continuously: 2=INC ONCE: ";B:IF B=0 OR B=1 THEN INPUT "INPUT T LIMIT,T1=";T1:GOTO 290
  30. 280 U=0:GOTO 350
  31. 290 IF (H>0 AND T1<T0+H) OR (H<0 AND T1>T0+H) THEN PRINT"T INPUT ERROR;":GOTO 270
  32. 300 GOTO 350
  33. 310 IF U=1 THEN 190
  34. 320 IF B=0 OR B=1 THEN GOSUB 600:GOTO 350
  35. 330 INPUT"MENU:0= CONTINUE:1= REPEAT last T increment :";U:IF U=0 THEN GOSUB 600
  36. 340 GOTO 190
  37. 350 GOSUB 370:GOTO 310
  38. 360 '----------------------------------------------MICROOTS to line 580
  39. 370 FOR J=0 TO N:R(J)=R:S1(J)=1:P(J)=0:Y(J)=C(J,0):NEXT J:'------  Y equivalent to X2,MICROOTS starting value/s
  40. 380 IF B=1 OR B=2  THEN PRINT TAB(3);"X"; TAB(25) "R";TAB(40) "P";TAB(50);"STEPS"
  41. 390 T=T0+H/2:FOR K=1 TO ST:L=-1:IF B=1 OR B=2 THEN PRINT TAB(49);K
  42.  
  43. 400 FOR I=0 TO N:GOSUB 680
  44. 410 ON A(I) GOSUB 10,11,12,13,14,15,16,17,18,19,20,21
  45. 420 IF ABS(Y(I))>1000000! OR ABS(F)>1D+30 THEN BEEP:PRINT"OVERLIMIT,see line 420.Repeat this T increment with different conditions:":K=ST:B=2:U=1:RETURN
  46. 430 Q=S(I)*ATN(F):S=SGN(Q)
  47. 440 IF S*S1(I)>=0 THEN P(I)=P(I)+G:R(I)=R(I)-G
  48. 450 R(I)=R(I)+G:Y(I)=Y(I)+Q*2^(P(I)/3-R(I)-1/3):IF B=0 THEN 470
  49. 460 PRINT Y(I);TAB(23);R(I);TAB(38);P(I)
  50. 470 IF ABS(Q*2^(P(I)/3-R(I)))<AC THEN L=L+1
  51. 475 A$=INKEY$:IF A$="P" THEN FOR J1=0 TO N:P(J1)=0:NEXT J1
  52. 480 A$=INKEY$:IF A$=" " THEN K=ST:B=2:U=1:RETURN
  53. 490 S1(I)=S:NEXT I:IF L<>N THEN 570
  54. 500 IF (H>0 AND T1<=T0+H) OR (H<0 AND T1>=T0+H) THEN B=2:' OR other ending conditions eg.  IF X(1,0)>*X(0,2)>0  THEN B=2
  55. 510 IF B=0 THEN BEEP:PRINT"T=";T0+H:K=ST:RETURN
  56. 520 BEEP:PRINT:PRINT"SOLUTION AT T=";T0+H:FOR I=0 TO N:PRINT"NEW C";I;"=";Y(I);
  57. 530 IF L=N THEN ON I+1 GOSUB 10,11,12,13,14,15,16,17,18,19,20,21
  58. 540 PRINT", RESIDUAL F";I;"=";F:NEXT I:PRINT K;" STEPS":IF B=1 THEN 560
  59. 550 IF E>0 THEN INPUT "REVIEW HIGHER ORDERS [eg.slopes X(I,1)], (1/0)=";K1:IF K1=1 THEN GOSUB 630
  60. 560 K=ST:RETURN
  61. 570 NEXT K
  62. 580 BEEP:PRINT ST;" STEPS,NO SOLUTION.Repeat this T increment with different conditions:":B=2:U=1:RETURN
  63. 590 '------------------------------UPDATE INITIAL CONDITIONS,to line 610
  64. 600 FOR I=0 TO N:C(I,0)=Y(I)
  65. 610 FOR J=1 TO O(I)-1:C(I,J)=2*X(I,J)-C(I,J):NEXT J:NEXT I:T0=T0+H:RETURN
  66. 620 '-------------------------REVIEW NEXT INITIAL CONDITIONS,to line 660 or functions thereof[ eg. X(0,0)*X(1,2)^3 ]
  67. 630 FOR I=0 TO N:PRINT"NEW C(";I;",";0;") will=";Y(I)
  68. 640 FOR J=1 TO O(I)-1:PRINT"NEW C(";I;",";J;") will=";2*X(I,J)-C(I,J):NEXT J
  69. 650 IF O(I)>0 THEN PRINT"NEW C(";I;",";O(I);") will=";X(I,O(I));"(approx,but not used)"
  70. 660 NEXT I:INPUT "PRESS RETURN to continue";C$:RETURN
  71. 670 '----------------------------DIFFERENTIAL SUBSTITUTIONS,to line 720
  72. 680 FOR A=0 TO N:IF O(A)=0 THEN X(A,0)=Y(A):GOTO 720
  73. 690 X(A,0)=(Y(A)+C(A,0))/2
  74. 700 FOR J=1 TO O(A):X(A,J)=2*(X(A,J-1)-C(A,J-1))/H:NEXT J
  75. 710 IF O(A)>1 THEN X(A,0)=(Y(A)+3*C(A,0)+H*C(A,1))/4
  76. 720 NEXT A:RETURN
  77. 730 '---VARIABLES LIMITERS when required for solution stability,to line 760.Tailor to the problem.
  78. '740 IF X(11,0)=0 THEN X(11,0)=.1
  79. '760 RETURN
  80.               MACINTOSH LISTING
  81.  
  82.